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Abstract 

We present a numerical evaluation of the loop-after-loop contribution to the 
second-order self-energy for the ground state of hydrogenlike atoms with low 
nuclear charge numbers Z. The calculation is carried out in the Fried- Yennie 
gauge and without an expansion in Za. Our calculation confirms the results 
of Mallampalli and Sapirstein and disagrees with the calculation by Goidenko 
and coworkers. A discrepancy between different calculations is investigated. 
An accurate fitting of the numerical results provides a detailed comparison 
with analytic calculations based on an expansion in the parameter Za. We 
confirm the analytic results of order a^(Za)^ but disagree with Karshenboim's 
calculation of the a^(Za)® In^(Za)^^ contribution. 



Typeset using REVTeX 



*e-mail: yerokhin@fn.csa.ru 



1 



INTRODUCTION 



In the \ow-Z region, calculations of radiative corrections in bound-state QED have his- 
torically relied on a (semi-) analytic expansion in powers of the external binding field Za. 
Calculations based on this perturbative approach have made an enormous advance during 
the last 50 years and achieved an excellent agreement with experiments (see, for example, a 
recent review [^). However, calculations in higher orders in Za become increasingly com- 
plex, as the number of terms in each higher order increases rapidly. Beside this, it is difficult 
to estimate the contribution of unevaluated higher-order terms. These are the reasons why 
the exact numerical treatment of radiative corrections is highly appreciated even in the low- 
Z region. It allows to test the reliability of methods based on an expansion in Za and can 
provide even more accurate results than analytic perturbative calculations. Some examples 
of this are the calculation of the self-energy correction to the hyperfine splitting in muonium 
performed by Blundell and coworkers 0, the calculation of the relativistic recoil correction 
for hydrogen by Shabaev et al. 0, and the evaluation of the first-order self-energy correction 
for Z = 1 — 5 by Jentschura et al. [|]. 

The aim of the present work is a numerical evaluation of the loop-after-loop contribution 
to the second-order Lamb shift of the ground state in hydrogen-like atoms to all orders in 
Za in the low-Z region. Analytic calculations of the Za-expansion coefficients for this con- 
tribution were previously performed by Eides and coworkers |^ and Pachucki |^ in order 
a^{ZaY and by Karshenboim in order a^{ZaY\r?{Za)~'^ . The first calculation of the 
loop-after-loop correction without an expansion in Za was carried out by Mitrushenkov et 
al. [§] for high-Z atoms. Recently, this correction was calculated to all orders in Za for the 
entire range of nuclear charge numbers by Mallampalli and Sapirstein 0. A fit to the data 
from Ref. confirms the analytic result of order a^{Za)^ but it is in a significant disagree- 
ment with Karshenboim's result of order a^ {Z aY \v? [Z a)~'^ . The subsequent calculation by 
Goidenko et al. [0, also non-perturbative in Za, shows to be compatible with the analytic 
calculations. In this work, we perform an independent calculation of the loop-after-loop cor- 
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rection and investigate possible reasons for the discrepancy between different calculations. 
Relativistic units are used in this article {h — c — m — 1). 



I. BASIC FORMALISM AND NUMERICAL PROCEDURE 

The expression for the irreducible contribution of Fig la (we refer to it as the loop-after- 
loop correction) reads 

where denotes the renormalized self-energy operator, \a) indicates the initial state and 
the summation is performed over the spectrum of the Dirac equation. The term with 
£„ = £„ corresponds to the reducible contribution and should be calculated together with 
the remaining diagrams in Fig. 1 . The self-energy operator is defined by its matrix elements 

/oo r 
du! / d^yiid^X2ipl{:s.i) 
-oo J 

xa^G{e - uj, Xi, X2)Q;^V6(x2)-D''''(a;, X12) 

-5m J c^3x^t(x)/5^,(x) , (2) 

where = (l,a); /3, a are the Dirac matrices, G{uj) = — — iO)) is the Dirac- 
Coulomb Green function, 7i = {cx -p) -\- j3m -\- Viii) is the Dirac-Coulomb Hamiltonian, 5m 
is the mass counterterm, and D^^(a;,Xi2) is the photon propagator in a general covariant 
gauge 



(3) 

kO=a; 



(27r)3 V k2 + i0 ^ '(k2 + i0)2 

To our knowledge, up to now all the practical self-energy calculations without an expansion 
in Za were carried out in the Feynman gauge (A = 1) which is technically the easiest 
choice of the gauge. While the usage of the Feynman gauge in calculations of the self- 
energy matrix elements is natural in the high-Z region, for low Z it is known to provide a 
spurious contribution of order Za which should be cancelled numerically to give a residual 



of order (Za)^. This spurious term is known to vanish in the Fried- Yennie gauge |1TT|JT2 



(A = 3), which possesses remarkable infrared properties. Since the present work is aimed 
to a calculation of the loop-after-loop correction in the low-Z region, we use the fact that 
this contribution is invariant in any covariant gauge and perform our calculations in the 
Fried- Yennie gauge. 

A general method which was used here for the calculation of the self-energy matrix 



elements can be found in Ref . , with some modifications due to a non-diagonal nature of 
the matrix elements and the different gauge. The self-energy matrix element is considered as 
a sum of two contributions originating from an expansion of the bound electron propagator 
in terms of interactions with the external field of the nucleus 

(a|S^(£)|6) = (a|4°+^)(5)|6) + {a\^^'^\s)\b) . (4) 

Here, the first term contains zero and one Coulomb interaction with the nucleus, and the 
second term contains two and more interactions. They are calculated in momentum and 
coordinate space, respectively. 

The expression (|I]) for the loop-after-loop contribution contains a summation of non- 
diagonal self-energy matrix elements over the whole spectrum of the Dirac equation. To 
perform the summation, we use the B-splines method for the Dirac equation developed by 
Johnson et al. [M . In this method, the infinite summation in the spectral representation of 



the Green function with a fixed angular momentum quantum number is replaced by a finite 
sum over basis-set functions. A straightforward evaluation of the sum in Eq. @j impl les 
a computation of many self-energy matrix elements with highly-oscillating wave functions 
and is computationally intensive. To reduce the computational time significantly, we define 
a self-energy correction to the wave function, as proposed in Ref. 

\^SE)=^R{ea)\a) . (5) 

According to Eqs. (H) and (|^), we write Eq. (|I]) as 

A TTi f '^^Pl '^^P2 (0+1) t/ N^rcd/ N (0+1)/ X 



+2 / ^ci3x,4°;^)'(Pl)G-'^(^a,Pl,X,)v.(¥(x2) 

+ J d'^,d'^2^^^f{^^)G'^\ea,^u^2)VsPi^2) , (6) 

where G"''^'^(£:a, xi, X2), G"''^'^(ea, Pi, P2), and G"'*^'^(eaj Pi, X2) are the reducible Dirac- Coulomb 
Green functions in coordinate, momentum, and mixed representations, respectively (by the 
mixed representation we mean the Fourier transform over one coordinate variable). 

As the first step of the numerical evaluation of Eq. (^, the effective wave functions 
(p^SE^\p) ^^'^ V^SE^x) are calculated on a grid and stored in an external file. Their com- 
putation is not much more intensive than an evaluation of a single self-energy matrix ele- 
ment. The most difficult part of the calculation is the evaluation of (p^g^\x.). Working in 
the Fried- Yennie gauge, we do not encounter severe cancellations between zero-, one-, and 
many-potential terms, as occur in the case of the Feynman gauge. Still, significant cancel- 
lations arise in the computation of the Green function G^^^-* which contains two and more 
interactions with the external field. In our implementation it is evaluated by a point-by- 
point subtraction of the two first terms of the Taylor expansion from the Dirac-Coulomb 
Green function (see Ref. for details) 



G^'^^\e,xi,X2) = G{e,xi,X2) - G{e,xi,X2)\z=Q - Z ^-^G{e,xi,X2) ^ 



(7) 



To control the cancellations which arise in the low-Z region, we monitor the corresponding 
Wronskian difference A^^^^^e) which can be calculated analytically {Af^(e) is the Wronskian 
of the solutions of the radial Dirac equation). Another numerical problem is the partial wave 
expansion. Its convergence is somewhat slower in the case of the Fried- Yennie gauge than in 
the Feynman gauge. In actual calculations we extended the summation up to sixty partial 
waves. It was performed before all numerical integrations were carried out. The remainder 
after the truncation of the sum was estimated taking into account the asymptotic behaviour 
of the expansion terms. Several checks were made of calculations of ^^se^\p) ^^se\^)- 
In one, we compared the diagonal self-energy matrix elements to the known results for the 
first-order self-energy contribution We also calculated the irreducible contribution 



to the self-energy correction to the hyperfine sphtting in H-hke atoms and found a good 
agreement with Ref. |0. 

In the next step, we perform the radial integrations in Eq. (|^). The Dirac-Coulomb 
Green function in coordinate space is evaluated using a finite basis set constructed from 
B-splines, after a transformation to a piecewise-polynomial representation as described in 
the Appendix. The momentum and the mixed representations of the Green function are 
obtained by the direct numerical Fourier transformation of the polynomial basis. After 
that, two-dimensional radial integrals in Eq. are expressed as a linear combination of 
one-dimensional integrals and can be easily evaluated up to a desirable precision. In actual 
calculations we used a basis set consisting of 70 positive and 70 negative energy states. The 
stability of the final results with respect to the size of the cavity and the number of energy 
states was checked. 



II. NUMERICAL RESULTS AND DISCUSSION 

In Table I and Fig. 2 we present the results of our calculation of the loop-after-loop 
contribution to the second-order Lamb shift of the ground state of hydrogenlike atoms, 
expressed in the standard form 

Ai?i,i= f-)'^Giai(Za) . (8) 

The results of two previous non-perturbative calculations of this correction are presented 
in Table I and Fig. 2 as well. A comparison exhibits a good agreement of the present 
calculation with the results of Mallampalli and Sapirstein and a strong deviation from 
the results of Goidenko et al. [jlO| . 



Let us consider possible reasons for this discrepancy. The method used in Ref. [|1^] is 
based on the multiple commutator approach combined with the partial- wave renormalization 
(PWR) procedure. In the PWR method, the truncation of the partial-wave expansion fulfils 
the role of the regularization parameter. This shows that this method is non-covariant. Still, 
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it can be used for the calculation of the diagonal first-order self-energy matrix elements, for 



which the PWR procedure is known to provide the correct result |l6| , p^ . In Ref. |T8[ this 
renormalization procedure was investigated for the self-energy correction to an additional 
Coulomb screening potential. It was shown analytically that some spurious terms arise in 
different parts of the total self-energy correction due to the non-covariant nature of the 
renormalization procedure. According to Ref. the spurious terms cancel each other if 



the perturbation is the Coulomb potential. The cancellation of the spurious contributions 
in the total self-energy correction holds no longer if the perturbation contains a magnetic 



photon (see Ref. [|T9[ and a conclusion remark in Ref. [T^ 

To consider this topic in more detail, we calculate the self-energy correction in the pres- 
ence of the perturbing potential —a/r both in the PWR scheme and using a covariant 
renormalization. For this choice of a perturbing potential, the total self-energy correction 
for a state \a) is d/ {dZ){a\T,ii{ea)\o) ■ The results of calculations are listed in Table II. Our 
calculation confirms the conclusions from Ref. |T^ about a) the presence of spurious terms 
in different parts of the correction and h) their cancellation in the sum for this particular 
choice of a perturbing potential. Summarizing, we conclude that it is possible that the 
PWR method applied to the irreducible part of the second-order self-energy correction, can 
provide a nonzero spurious contribution. 

In order to compare our results with calculations based on an expansion in Za, we 
approximate our data for the function Giai by a least-squares fit with five parameters 050, 
o-mi '^625 ciQii and ago (the first index of the a coefficients indicates the power of Za, the 
second corresponds to the power of ln(Za)~^). A fit to our numerical results in Table I 
yields 

aso = 2.33 aes = -1.1 ■ (9) 

This is in a good agreement with the fitting coefficients from Ref. [^] (050 = 2.3 or 2.8 for 
different sets of data, = —0.9) but disagrees significantly with Karshenboim's analytic 
result = -8/27 0. 



In order to investigate this discrepancy in more detail, we note that the Za-expansion 
calculations of the loop-after-loop correction in Refs. 1^-0 were performed in the Fried- 
Yennie gauge like in the present work and, therefore, it is possible to compare the calculations 
on intermediate stages. So, we expand the inner electron propagators in diagram Fig. la 
in terms of interactions with the nuclear binding potential and calculate the first six terms 
of the expansion separately. The corresponding Feynman diagrams are presented in Fig. 3. 

(2+) 

These diagrams do not contain y^^^ (x) which is the most difficult part of the calculation. 
Therefore, we were able to calculate them for very low fractional Z . This is important for 
a reliable fitting of our data which vary very fast in the vicinity of Z = 0. The remainder 
behaves more smoothly in the low-Z region and its fitting is easier. In the calculation of 
the diagrams shown in Fig. 3, we use closed analytical expressions for the Dirac Green 
function with zero and one Coulomb interaction. In this way we eliminate the numerical 
uncertainty due to the finite basis set representation of the Green function. The numerical 
results for each diagram in Fig. 3 were approximated by least-squares fits with eight or 
seven parameters 050, agj (i = 3, . . . , 0), a^i {i = 3, 2, 1) (in the last case 071 was omitted). In 
order to reduce the statistical uncertainty of the fitting procedure, a large number of points 
(twenty or more) was used. The stability of the fitting coefficients was checked with respect 
to the number of points, minimal and maximal nuclear charge numbers, and different fits. 
The numerical results and the fitting coefficients for diagrams in Fig. 3 are listed in Table 
III. 

We found a good agreement with results from Refs. for the coefficient 050 and with 
Ref. for the coefficient ags originating from diagram Fig. 3f. The only discrepancy with 
the analytical calculations originates from diagram Fig. 3c. While this diagram should not 
contribute to order a'^{Zay \n'^(Za)^'^ according to Karshenboim, our calculation shows the 
presence of a cubed logarithm with coefficient aga = —0.652(30). 

Summarizing, we conclude that our calculation of the loop-after-loop correction confirms 
the analytic result of Refs. [^,0 for the coefficient 050 (050 = 2.3). A fit to the numerical 
results yields 
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= -0.958(30) 



= 3.3(5) 



(10) 



for the diagrams shown in Fig. 3, and 



aes = -0.05(7) 



a62 = 1.2(8) 



(11) 



for the non-perturbative remainder. 

We note a remarkably slow convergence of the Za-expansion for the loop-after-loop 
contribution to the second-order Lamb shift. As an illustration, in Fig. 4 we plot the contri- 
butions of the first one, two, and three expansion terms together with the non-perturbative 
results. The expansion coefficients are taken from Eqs. (0) and (pT|). One can see that even 
for hydrogen the contribution of the first three expansion terms covers only about 50% of 
the total result. To obtain a reasonable fit to the numerical data even for very low Z, it 
is necessary to take into account at least four first expansion terms. This fact shows the 
necessity for non-perturbative (in Za) calculations of the total second-order Lamb shift in 
the low-Z region. 
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TABLES 

TABLE L The loop-after-loop contribution to the second-order Lamb shift of the ground state 
of hydrogenlike atoms expressed in terms of the function Giai(.Z'a) defined by Eq. @. 
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-2.694 


8 


-4.998(1) 




-2.659 


9 


-4.958(1) 




-2.642 


10 


-4.902(1) 


-4.9016(14) 


-2.601 


12 


-4.762(1) 








15 


-4.523(1) 


-4.5218(6) 






20 


-4.122(1) 


-4.1217(3) 


-2.568 
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TABLE II. The self-energy correction in the presence of the perturbing potential — a/r cal- 
culated both in the partial-wave renormalization scheme (PWR) and using a covariant renormal- 
ization (CR). A£?ir is the irreducible contribution (known also as perturbed orbital contribution), 
A^vr denotes the sum of the reducible and the vertex contribution. The results are compared with 
the derivative of the first-order self-energy contribution AEse with respect to the nuclear charge 
number Z. The calculation is performed in the Feynman gauge for a point nucleus. 



z 




AEP^ 

ir 


A^PWR 








AEsE/dZ 


20 


0.02003 


0.00813 


-0.00899 


0.00289 


0.01104 


0.01102 


0.01102 


30 


0.03969 


0.02137 


-0.01080 


0.00750 


0.02888 


0.02886 


0.02885 


50 


0.10722 


0.07372 


-0.00892 


0.02460 


0.09830 


0.09832 


0.09832 


70 


0.23749 


0.18248 


0.00197 


0.05698 


0.23946 


0.23946 


0.23946 


92 


0.55983 


0.46257 


0.03473 


0.13199 


0.59456 


0.59456 


0.59456 
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TABLE III. The contributions of the diagrams shown in Fig. 3, expressed in terms of the 
function G\i^i{Za) defined by Eq. (|8|). The numerical results for the first three coefficients of 
the Za-expansion corresponding to two different fits are listed and compared to the analytical 
calculations. 



z 


Fig. 3a 


Fig. 3b 


Fig. 3c 


Fig. 3d 


Fig. 3e 


Fig. 3f 


Fig. 3(a-f) 


0.1 


0.0072 


9.1740 


-7.9975 


-0.0099 


0.1654 


-0.7931 


0.5460 


0.2 


0.0121 


9.0885 


-8.4197 


-0.0176 


0.2709 


-1.1941 


-0.2599 


0.4 


0.0198 


8.9417 


-8.9281 


-0.0306 


0.4334 


-1.7448 


-1.3086 


0.7 


0.0283 


8.7535 


-9.3648 


-0.0472 


0.6204 


-2.3078 


-2.3176 


1.0 


0.0347 


8.5886 


-9.6271 


-0.0618 


0.7705 


-2.7177 


-3.0128 


1.5 


0.0425 


8.3471 


-9.8763 


-0.0831 


0.9730 


-3.2199 


-3.8168 


2.0 


0.0478 


8.1349 


-9.9989 


-0.1019 


1.1372 


-3.5886 


-4.3696 


3.0 


0.0534 


7.7705 


-10.0572 


-0.1346 


1.3949 


-4.1007 


-5.0738 


5.0 


0.0526 


7.1957 


-9.8680 


-0.1876 


1.7519 


-4.6721 


-5.7274 


7.0 


0.0425 


6.7505 


-9.5453 


-0.2309 


1.9946 


-4.9565 


-5.9451 


10.0 


0.0177 


6.2353 


-9.0170 


-0.2852 


2.2458 


-5.1386 


-5.9419 


15.0 


-0.0376 


5.6313 


-8.2016 


-0.3600 


2.5111 


-5.1632 


-5.6201 


20.0 


-0.1023 


5.2261 


-7.5211 


-0.4248 


2.6849 


-5.0591 


-5.1963 


Analytic 


results 














^50 





9.284 


-6.984 











2.300 


am 

















-0.296 


-0.296 


Eight-parameter fit: 














a50 


0.000 


9.284 


-6.985 


0.000 


0.000 


0.001 


2.300 


063 


0.000 


0.001 


-0.658 


0.002 


-0.003 


-0.304 


-0.963 




0.02 


-0.09 


3.1 


-0.07 


1.17 


-0.75 


3.34 


Seven-parameter fit: 














^50 


0.000 


9.285 


-6.987 


0.000 


0.000 


0.000 


2.298 
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062 



-0.001 
-0.01 



-0.003 
-0.01 



-0.646 
2.82 



0.002 
-0.08 



-0.005 
1.22 



-0.301 
-0.81 



-0.952 
3.12 



APPENDIX: PIECEWISE-POLYNOMIAL REPRESENTATION OF THE GREEN 

FUNCTION 

The B-splines method for the Dirac equation provides a finite set of radial wave 



functions with a fixed angular momentum quantum number which can be written in the 
form 

^Ux) = -Y.c',{^,n,l){x-xi)' , (Al) 

^ k 

assuming that x G [x/, x^+i]. Here xi is the radial grid, the index i = 1,2 indicates the upper 
and the lower component of the radial wave function, n numbers the wave functions in the 
set, and k is the angular momentum quantum number. The radial Green function, defined 
by 

G«(.,.„.,) = i:<^i(iiM^, (A2) 
can be written in the piecewise-polynomial representation as follows: 

G'^{e,x,,X2) = ^ E Al,^{e,K,h,k){x-xi,)'^{x-xi,)'^ , (A3) 
where Xi G [xi-^,xi-^+i], X2 G [xi.^,xi.^^i]. The coefficients ^^"'^j^^ 

The radial Green function in momentum space can be written in the same way using 
Fourier transformed basic polynomials 

Gne.p1.p2) = E E Kk.{e,KM,hwt:{pi)^i'{P2) , (A5) 

nf (p) = AnsiLi) f dxx{x - xifjL^P^) , (A6) 

Jxi 

where Li^2 = |k ± 1/2| — 1/2; s{Li) = 1, 5(^2) = and jiiz) denotes the spherical 

Bessel function. 
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FIGURES 





a b 
FIG. 1. One-electron self-energy Feynman diagrams of second order 
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FIG. 2. The function G\a,i{Za) in different calculations. The solid line indicates a fit to our 
numerical results. 
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=c^^o= 



a b(2) 



I II III 

^ <C> <0> <C> ^ <5> 



d e(2) f 

FIG. 3. Diagrams obtained from Fig. la by expansion of the inner electron propagators in 
terms of interactions with the nuclear binding potential. A double line denotes the electron in the 
field of the nucleus. A single line indicates the free electron. A dashed line denotes a Coulomb 
interaction with the nucleus. Some diagrams are counted twice, as is denoted by "(2)". 
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FIG. 4. The non-perturbative (in Za) function G\s\{^Za) and the contributions of the first one, 
two, and three terms of its expansion in Za. Dots indicate the non-perturbative results. The 
expansion coefficients are taken from Eqs. ( p!o| ) and (11). 
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